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ABSTRACT 


-J 

A  general  theory  for  nonlinear  implicit  one-step  schemes  for  solving  initial  value  problems  for 
ordinary  differential  equations  is  presented  in  this  paper.  The  general  expansion  of  ’symmetric’ 
implicit  one-step  schemes  having  second-order  is  derived  and  stability  and  convergence  are  studied. 
As  examples,  some  geometric  schemes  are  given. 

Based  on  previous  work  of  the  first  author  on  a  Generalisation  of  Means,  a  fourth-order 
nonlinear  implicit  one-step  scheme  (GMS)  is  presented  for  solving  equations  with  steep  gradients. 
Also,  a  hybrid  method  based  on  the  GMS  and  a  fourth-order  linear  scheme  is  discussed.  Some 
numerical  results  are  given. 


1.  Introduction 


Many  classical  methods  for  solving  initial  value  problems  for  ordinary  differential  equations  are 
based  on  piecewise  polynomial  interpolation.  If  the  solution  of  the  problem  possesses  a  very  steep 
gradient,  these  schemes  produce  poor  results.  In  particular,  if  a  singularity  occurs,  it  is  often 
inappropriate  to  attempt  to  represent  the  solution  in  the  neighborhood  of  the  singularity  by  a 
polynomial.  In  this  paper,  we  consider  a  class  of  nonlinear  implicit  one-step  schemes  that  may  be 
more  appropriate  for  such  problems. 

A  genera]  theory  for  nonlinear  implicit  one-step  schemes  is  developed  in  section  2.  Conditions 
for  consistency,  stability,  and  convergence  are  obtained.  Each  consistent  symmetric  scheme  is  at 
least  second-order,  and  the  condition  that  it  must  satisfy  to  be  fourth-order  is  given.  A  class  of 
symmetric  and  homogeneous  schemes  which  are  generalisations  of  tEe  well-known  trapezoidal  rule  is 
obtained. 

The  trapezoidal  rule  is  exact  for  second-degree  polynomials.  In  terms  of  geometry,  a  second- 
degree  polynomial  is  a  conic.  As  examples  of  nonlinear  symmetric  implicit  schemes,  we  develop 
several  geometric  schemes  based  upon  ’’circles",  "ellipses",  "parabolae",  and  "hyperbolae”  in  section 

3. 

On  the  other  hand,  in  terms  of  Means,  the  trapezoidal  rule  is  the  Arithmetic  Mean  of  the  first 
derivative  of  the  solution  at  two  neighboring  grid  points.  In  section  4,  based  on  the  Generalization  of 
Means  [9],  a  fourth-order  nonlinear  implicit  one-step  scheme  (GMS)  is  presented  for  solving  problems 
with  steep  gradients. 

In  section  5,  we  discuss  some  practical  considerations  including  the  use  of  hybrid  methods  based 
upon  the  GMS  and  more  traditional  schemes. 
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Id  this  paper,  the  theory  of  nonlinear  implicit  one-step  schemes  is  restricted  to  scalar  equations. 
However,  we  have  used  these  schemes  successfully  to  solve  systems  of  equations.  The  application  of 
these  schemes  to  systems  is  discussed  briefly  in  section  6. 

Numerical  Results  for  seven  test  problems,  some  of  which  contain  systems  of  equations,  are 
given  in  the  last  section.  Two  of  the  examples  use  an  imbedding  technique  to  apply  the  GMS  to  the 
solut  ion  of  two-point  boundary  value  problems. 
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2.  A  general  theory  for  nonlinear  implicit  one-step  schemes 

Consider  the  initial  value  problem  (I.V.P.) 

y’ —  f(x,y),  y(a)  — ya  (a<x^.b)  (1) 

where  f(x,y)  is  continuous  in  x  and  Lipsehitz  continuous  in  y  in  the  region 
a^x<£b,  -  oe  <  y  <  oo,  a  and  b  finite. 

We  investigate  the  following  general  nonlinear  implicit  one-step  scheme 

Yn+1  “  Yn  +  b  S(fn,fn+1),  (2) 

where  ' 

**  *  Xn+1  ’  Xn’  “  ^Xn’^  ^n+1  ™  ^Xn+l’^n+l^ 

The  the  local  truncation  error  for  scheme  (2)  is 

m  __ 

L(f)  «  y(xn+,)  -  y(xn)  -  hS(f(xn,y(xn)),f(xB+1,y(xii+1))),  (3) 

where  y(x)  is  the  solution  of  (1). 

Definition  1:  [o]  The  scheme  (2)  is  said  to  have  order  p  if  p  is  the  largest  integer  for  which 

L(f)  -  0(hP+1). 

Definition  2:  The  scheme  (2)  is  said  to  be  consistent  with  the  I.V.P.(1)  if  L(f)  ■«*  o(h). 

We  will  use  the  notation  f(t)mf(t,y(t))  throughout  this  paper  except  where  it  may  be  confused. 
For  xn  £.x  ^Lxn+1,  let  t  »■  (x  -  xn)/h.  Since 
y(xn+i)  *  y(xn)  **  h/J  f|t)  dt, 

(3)  may  be  rewritten  as 

Uf)  —  h{/J  f(t)  dt  -  S(f(xn,y(xn)),f(xn+1,y(xn+1)))}.  (4) 

By  the  Integral  Mean  Value  Theorem,  there  exists  a  point  (  between  xn  and  xn+ ,  such  that 

/  ’  f(t)  dt-ftfjiO). 

o 


So 


tP-  !tf,y(0)  •  S(f(xn,y(xn)),nxn+1,y(xn+1))). 
df 

Furthermore,  if  T(t)  —  exists,  then 

/‘imdt-  nxj  +  z'n txi-t)dt 

and 

S(f(xn),f(xn+1))  -  S(f(xn),ffxn))  +  /o‘  dS(f(xn),f(t)). 

Hence,  we  have 

Lemma  3:  Let  the  function  S(f,g)  be  continuous  in  its  two  variables  f  and  g.  Then  the 
scheme  (2)  is  consistent  with  (1)  if  and  only  if 

S(f,f)  -  f.  (5) 

Moreover,  if  both  f  and  S  have  continuous  first  derivatives,  then  any  scheme  (2) 
satis  fying  (5)  has  truncation  error 

W)  -  h/o‘{r(t)(l-t)  -  S’(f(xn,y(xn)),f(t))}dt,  (6) 

where  f’(t)  -  ^  S’(g.flt))  -  jM! 

The  proof  of  the  following  theorem  is  similar  to  the  one  given  in  [2]  for  general  explicit  one-step 
methods. 

Theorem  4:  Let 

(i)  the  function  S(f,g)  be  continuous  jointly  as  a  function  of  its  two  arguments  in  the 
region  fg  ^  0,  and 

(ii)  S(f(z,y),g(x,y))  satisfy  the  Lipschitz  condition 

|  S(f(x,y),f(i,w))  -  S(f(x,y*),f(i,w*))|  ^  M  (  |y  -  y*|  +  |w  -  w*|  ) 

for  all  points  in  the  domain  defined  by 

>  *  * 
a  <£.  x,t  <L  b,  -oo  <  y,y  ,w,w  <  +oo 

under  the  constraints  f(x,y)f(z,w)^p  and  f(x,y*)f(z,w*)^>0. 

Then  the  scheme  (2)  is  convergent  if  and  only  if  it  is  consistent. 


In  order  to  get  a  second-order  scheme,  note  that 
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/o1  f(t)  at « nxn)  +  |f’(xn)  +  ntxi-t)2  at, 

lol  dS(f(xn),f(t))  -  h  |jS(f(xn),f(x))|x_x>  +  /0‘  S”(f(xn),f(t))  (l-t)dt. 

Therefore, 

Theorem  5:  If  S,f  t  C1,  then  the  eeheme  (2)  has  a  second-order  rate  of  convergence  if 
and  only  if 


S(f.f) 


t  dS(g,f) 

'■  THr., 


1 

T 


(7) 


Moreover,  i f  the  second  derivatives  of  S  and  f  are  continuous,  then 


L(f)  -  -  ][  /  '{(l-t)2r(t)  -  2(l-t)S”(f(xn),f(t))>  dt 


where  S(f(zJ,f(t))  —  S( f(xn,y( zj),f(t,y(t))). 

Corollary  6:  Keeping  S(f,f)"*f,  the  second  condition  of  (7)  is  equivalent  to  one  of  the 
following  four  conditions 


dS(g.f)  _  1  3S(g,f)  _  1 

iri,.,  r  iflH “r 


gsig.n 


as(g.t) 


■ar"W~srw  “d  TT^T 

Suppose  now  S  and  f  both  belong  to  C2.  Expanding  f  about  the  point  \+1/2  “  (xn+xn+,)/2 
we  have 

r  fit)  dt  -  fn+1/! + ^r’(x„tl/!) + o(h«). 

On  the  other  hand, 

dSlflRfln)) 

(Vxn^l/2)  + 


gS(g,f) 


(8) 


5S(f(e),f(«7)) 

S(Rxn),f(xn+1))  —  S(^xn+l/2^xn+l/2W  + 


dq 


a’sfflOfW)  (Vx.+I  It?. 

dft" — '  ■  ■  ■  ■  ■ 


(xn+l*xn+l/2^  +  qM  I 


i-fl/2 


dHmm) 
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(xn+l-xn+I/2KVx„+l/2)  +  —To 


d2mmn)) 


-+ 


For  symmetric  schemes,  S(f,g)*»S(g,f)  and,  consequently, 

h2  ^S(f(0,f(9)) 

s(f(Xn),f(xn+,))  -  f(xn+1/2)  +  jp{  — - + 


d2S(mf(9))  d‘S(f{af(9))  tl  , 

s? — 2-^r1u^b) 


d2sm,t(n))  as(f(C),f|t>)L .. 


d*sm,f(v))  a2s(f(fl,f(g)) 

dt  dt,  "  WvMO 

Hence,  from  Corollary  6 


nma 


h2  .  d*S(f,g)  d*S(f,g), 

sM*.M(*.+1»  -  «*„♦„.)  +  it  r  +  .  _ _ ] 


)  +  0(h«). 


Substituting  (9)  and  (8)  into  (4),  we  obtain  the  following  theorem. 

Theorem  7  s  If  S  and  ft  Cs,  then  each  symmetric  consistent  scheme 

S(f,f)  -  f,  S(f,g)  -  S(g,f)  (10) 

has  a  second-order  rate  of  convergence,  at  least.  Moreover,  if  the  fourth  derivatives  of  f 
and  S  are  continuous,  then 

h3  o  c>2S(f,g)  d2S(f,g)  . 

utH .  ^  (rw^-jjs-  -  )+  o(bs,.  (») 

It  should  be  noted  that  symmetry  is  a  sufficient  but  not  necessary  condition  for  a  scheme  to  be 
second-order.  For  example,  the  scheme 

Y  Y  ■  , 

otl  "  '+V.+W/* 

is  second-order  but  not  symmetric. 

Corollary  8:  A  symmetric  scheme  is  fourth-order  if 


»   st*'.'  w  -  1 


a2S(f,g)  d2S(f,g) 


9f" 


3f3g 


T’ 

^  ^.S“f(*1+i/2)  3f2  ^X“*»+l/2 


(12) 


Now  we  consider  a  general  representation  of  consistent  symmetric  schemes  (2).  Let 
e  _  f*g  _  „  f+g 

and  assume  S(f,g)  can  be  expanded  in  terms  of  its  two  variables  f  and  g: 

S(f,g)  *  S(£,ij)  —  E  ok  j  {k 


k.i 


Using  conditions  (10),  each  consistent  symmetric  implicit  one-step  scheme  (2)  has  the  expansion 
S(f,g)  -  ^  +  E  a2kJ  (^)2k+2  (£^»,  (13) 


»2P.j 

where  a0.  ■  are  real  constants  to  be  chosen. 

*K,J 

Furthermore,  it  is  often  useful  to  restrict  the  class  of  schemes  to  be  homogeneous  in  the  sense 

that 

S(cf,  eg)  —  cS(f,g) 

for  any  constant  c.  For  these  schemes,  we  obtain  the  following  conclusion: 

Theorem  9:  Assume  S  can  be  expanded  in  terms  of  its  two  variables  f  and  g.  Then  each 
homogeneous  consistent  symmetric  nonlinear  implicit  one-step  scheme  (2)  has  the 
following  expansion 

f+g  ,  .  ^  _  /f-g  »2k+2 


(14) 


S(f,g)«^{l-  E  «k&2k+M- 

where  ak  are  real  constants  to  be  determined. 


(15) 


Observe  that  the  trapezoidal  rule  is  the  principal  part  of  each  homogeneous  consistent  symmetric 
nonlinear  implicit  one-step  scheme  (2).  Hence,  in  this  sense,  these  nonlinear  schemes  are  an  extension 
of  the  trapezoidal  rule. 


Setting  crk  ■=  0  for  all  k  >  0,  we  get  an  'extended  trapezoidal  rule’  with  one  extra  term: 
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From  (11), 


h3  r- 

**  f  O -■ 


L( f;  o)  =  -  ^  f”  -  30^}  +  0(ha).  (17) 

Id  terms  of  Means,  the  scheme  (16)  represents  a  linear  combination  between  the  Arithmetic 
Mean  and  the  Harmonic  Mean  of  f  and  g: 

S(f,g)  =  (l*a)  +  Q 

For  example,  if  a  =  1,  the  above  scheme  represents  the  Harmonic  Mean  between  f  and  g. 

Finally,  we  discuss  the  stability  of  the  nonlinear  implicit  one-step  scheme  (2). 

Let  f(x,y)  =  Xy  and  Yn+1  *■=  pYn,  where  p  is  the  "growth  factor"  in  the  step.  Assuming  that  the 
scheme  (2)  is  homogeneous  in  the  sense  of  (14),  we  get 

p  *  1  +  hX  S(l,  p).  (18) 


Definition  10:  A  nonlinear  implicit  one-step  scheme  (2)  is  said  to  be  A-stable  if  all  the 
roots  of  its  characteristic  equation  (18)  satisfy  \p\  <  1  for  any  Re  hX  <  0. 


For  nonlinear  implicit  one-step  schemes,  (18)  may  have  more  than  one  root  for  a  fixed  hX,  and  it 
may  be  possible  to  choose  which  root  the  scheme  follows,  unlike  the  case  for  multistep  methods. 

Hence,  the  following  definition  may  be  of  some  practical  value. 

Definition  11:  A  nonlinear  implicit  one-step  scheme  (2)  is  said  to  be  conditionally  A-stable 
if  at  least  one  root  of  its  characteristic  equation  (18)  satisfies  \p\  <  1  for  any  Re  hX  <  0. 

Theorem  12:  For  each  real  symmetric  homogeneous  scheme  S(f,g),  the  characteristic 

equation  (18)  transforms  the  unit  circle  of  the  p  plane  to  the  imaginary  axis  of  the  X 

plane. 

Proof:  Since  S(f,f)  —  S(f,f)  *  S(f,f),  S(f,f)  is  real  for  any  f.  Hence,  for  p  =  e1^, 
e*-l  e^-eW2 

Xh  «  — rr-  *  — rrys — is  purely  complex. 

S(e'^,  1)  S(e‘^2,  e"^2) 

A  necessary  requirement  for  a  homogeneous  scheme  to  be  A-stable  is  that  it  is  stable  at  infinity. 


Theorem  13:  A  necessary  condition  for  a  homogeneous  nonlinear  implicit  one-step 
scheme  (2)  to  be  A-stable  is  that  all  roots  of 

S(l,  />)=-=  0  (19) 

satisfy  |i|  1.  To  be  conditionally  A-stable,  at  least  one  root  of  the  above  equation  must 

satisfy  |i|  £  1. 

As  an  example,  consider  the  stability  of  scheme  (16)  with  characteristic  equation 

,  _  ,  +  bx  (£+1-  «£$ ).  (20) 

If  |i|  =  1,  then  i  =*  e1^,  0£j)<2n,  and 

JTl  -  i  tan  m- 

Hence,  (20)  may  be  rewritten  as 
2i  tan  0/2 

h\  * - t— 

l+o[tan  0/2]' 

It  follows  that,  if  |i|  *=  1,  then  hX  is  purely  complex  and  lies  in  the  interval  (-ia*1/2,  iof1/2)  for 
a  >0.  For  a<0,  hX  may  assume  any  value  on  the  imaginary  axis.  Also  note  that,  for  scheme  (16), 
equation  (19)  becomes 


That  both  roots  of  S(l,i)  are  on  the  unit  circle  for  a<0  and  that  one  is  inside  and  the  other  is 

outside  the  unit  circle  for  o> 0  follows  from  the  well-known  result 
Lemma  14:  The  one-to-one  mapping  in  the  complex  field 

W(z)~rin  (21) 

maps  the  domains  |s|<l,  |«|“/,  and  |«|>i  onto  ReW<0,  ReW—0,  and  ReW>0, 
respectively. 


The  characteristic  equation  (20)  can  be  rewritten  as 

®W2(i)  +  Lw(i)  -  1  -  0. 
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•where  W  is  defined  in  (21). 

Since  W=0  is  not  a  root  of  equation  (22),  any  root  W  satisfies  hX  {W1  +  (-o)W)  *  2.  If  o<0 
and  ReXh  <  0,  then  ReW<0,  whence,  by  Lemma  14,  any  root  of  the  characteristic  equation  (20)  for 
the  scheme  (16)  satisfies  \p\  <  1. 


Also  note  that,  if  oy^O,  then  the  roots  of  quadratic  equation  (22)  satisfy  WjW2  —  -  o’1.  Hence, 
if  o>0,  then  Arg(Wj)  +  Arg(W2)  —*  it  and,  consequently,  either  ReWj—ReWj—O  or  ReWj  and 
Re\Vr2  are  opposite  signs.  Therefore,  by  Lemma  14,  either  both  roots  of  the  characteristic  equation 
(20)  for  the  scheme  (16)  satisfy  |^|«-1  or  one  is  larger  than  1  in  magnitude  and  the  other  is  smaller. 
We  may  simply  choose  the  value  of  Yn+1  in  scheme  (2)  with  (16)  such  that  ||Yn+1||  <  ||Yn||.  Thus, 
we  have  proven 

Theorem  15:  Scheme  (16)  is  A-stable  for  a  £0  and  conditionally  A-etablc  for  a  >  0. 


For  the  more  general  scheme  (15)  with  a  finite  terms  number  of  terms,  the  corresponding 


characteristic  equation  is 

b  ok  (w(,)}2k+2  +  j^w(,)  -  i  -  o. 

By  the  relationship  between  coefficients  and  roots, 


E  W, 


-1  _2 


EX’ 


(23) 


(24) 


where  Wk  (k=l,...,2n+2)  are  roots  of  (23).  Hence,  if  Re(hX)  <  0,  then  ReWk < 0  for  at  least  one 
root  \Vk  of  (23).  Therefore, 

Theorem  10:  Each  scheme  (15)  with  a  finite  number  of  terms  is  conditionally  A-etablc. 


Remark:  Theorem  16  is  valid  even  if  the  coefficients  ok  of  the  scheme  (15)  are  complex. 
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3.  Derivation  of  some  geometric  schemes 

The  trapezoidal  rule  can  be  viewed  as  the  Arithmetic  Mean  of  fn  and  fn+l  since 

*«,«>  -  k + 

Let  fn  =  tan  an,  fn+,  —  tan  on+1,  (y(xn<H)-y(xn))/h  —  tan  an+1/r  The  trapezoidal  rule 
satisfies 

tan  <*n+1/2  —  j(tan  on  +  tan  an+1).  (25) 

It  is  easy  to  see  that  the  scheme  (25)  is  poor  if  the  angle  on  or  on+J  is  close  to  90°.  In  this  case, 
it  is  natural  to  replace  (25)  by  the  Arithmetic  Mean  of  the  angles  an  and  an+J: 

°n+l/2  “  +  an+l>-  (26) 

The  corresponding  function  S(fn,  f„+»)  is 

{(i+fn2xi+f„ .  ,2)}l/2  +  f  f. .  - 1 

at  t  \ _  '  n  n+l  "  n  n+I  _  io7» 


S(fn’fn+1^ 


f»  +  fn+I 


From  analytic  geometry,  the  curve  which  satisfies  (26)  everywhere  is  a  circle.  So  we  call  (27)  a 
Circle  Scheme.  The  Circle  Scheme  (27)  is  not  linear  with  respect  to  the  solution  y(x)  or  f,  but  (26)  is 
linear  with  respect  to  the  angles.  Hence,  if  the  angles  are  not  too  large,  the  Circle  Scheme  is  eloae  to 
being  linear.  In  fact,  if  we  rotate  the  coordinate  system  by  an  angle  P  =  an+,y2,  then  the  Circle 
Scheme  coincides  with  the  trapezoidal  rule  in  the  new  coordinates. 

Introducing  a  parameter  a  into  (27)  leads  to  a  class  of  Elliptic  Schemes: 

{(a2+f  2)(a2+f  2)}1/2  +  f  f  , ,  -  a2 

E( f  n’fn-*- 1  *a )  - Z - 7 - 7 - —  (28) 


W,g)  (a2+g2)‘/2  E  3E(f,g)_  (aM2)1/2  E 
df  (a2+f2)1/2(  f-t-g)’  dg  (a2+g2)1/2(f+g)’ 
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d-E  3EE-2f  d2E  dE  E 

df‘  df  a2+f2  dfdg  dt  a2+f2 


Cj  I  U  £j  1 

“  '  ^2^2,  ’  “  2(a2+f2)  ' 

A  straightforward  computation  leads  to  the  following  conclusion. 

Theorem  17:  E(fn,fn+J;a)  in  (28)  has  the  following  properties: 

(fn+fn+1)E^0  with ’-’iff  fn+fa+1-0. 

f-io.  *0.  (f.+f.+1)gio. 

n  n*f  I 

Min(fn  •W  E(r„.f„t,»)  £M«<f„,f„+l). 

if  .  >  (f„f„,),/!. 

'n^'n+l 

A*  a  function  of  a,  E  has  only  one  fixed  point  a  =  (fnfn+1t^>  for  fnfn+t>0. 

\ 

Hence,  the  Elliplic  Scheme  (28)  represents  a  Mean  which  lies  between  the  Arithmetic  Mean  and 
the  Harmonic  Mean. 


Similarly,  we  can  derive  two  other  geometric  schemes:  the  Parabolic  Scheme  and  the  Hyperbolic 
Scheme.  An  easy  way  to  derive  the  Parabolic  Scheme  is  to  apply  the  trapezoidal  rule  in  a  coordinate 
system  rotated  by  an  angle  a=*arctan(o)  from  the  original  coordinate  system: 


Substituting  hyperbolic  functions  into  the  formulas  (25), (26)  instead  of  trigonometric  functions,  we 

get  the  Hyperbolic  Scheme 

a2+f  f  . ,  -  {(a2-f n2)(a2-f„  2)},/2 

H(fn.fn+1:a) - ""-tL  -S+lZL  ,  (30) 


fn+fn+l  fn^n+l+a2 
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4.  Generalized  Mean  Scheme  (GMS) 

In  addition  to  the  above  Geometric  Means,  another  useful  Mean  for  solving  O.D.E.S  is  the 

Generalized  Mean  developed  by  Jiachang  Sun  (9]. 

Definition  19:  For  a  given  positive  sequence  a  “  (aj,...,an)  on  a  real  plane  (r,t),  a 
Generalized  Mean  of  the  sequence  {a}  S(aj,...,an;r,t)  is  defined  by 

(n-1)!  T(t+1) 

a. . .  -  (  ’  ,yw)V“,  (31) 


r  (t+n) 


where  [>1,...,yn]f(y)  is  the  (n-l)-th  divided  difference  of  the  function  f(y)  at  the  points 

Now,  we  use  the  Generalized  Mean  (GM)  in  (31)  to  construct  the  Generalised  Mean  Scheme 
(GMS),  a  nonlinear  implicit  scheme.  In  this  paper  we  only  consider  the  one-step  case.  From  (31), 
the  GM  between  fn  and  fn+1  is 


f  r(l+t)  .  f  r(l+t) 

- ( L  >'/n 


where  r,t  are  real.  Substituting  (32)  into  the  local  truncation  error  formula  (11),  we  get 

Theorem  20:  Let  f(x)  have  constant  sifn  for  Z(t  <  z  <  xn+J,  then  each  scheme  (32) 
with  two  real  parameters  (r,t)  is  second-order  at  least.  Moreover, 

Uf;  r)  -  -  £  {  r  -  [3  -  r(2+t)]^}|x  _  x>  +  h  +  0(h5),  where  h  -  xn+,  -  x„.  (33) 

To  simplify  the  study  of  this  scheme,  we  consider  the  restriction  rt=l  on  the  parameters  r,t. 

The  scheme  (32)  reduces  to 

f  1+r  _  f  1+r 

c it  r  _  t  n+l  n 


where 


S(f  f  Q)  —  n-fl  n  s(f  f  •-!)* _ n  n+ 1 

^fnJn+rU'  |  „,f  ^VSi+l’1'  f  -1  f -1’ 

'*n+l' ^n'  fn+l  *  fn 


and  the  local  truncation  error  (33)  becomes 


Uf;r)»  -^{r-d-r)^}  +0(h5). 


(36) 
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it  i$  obvious  that  for  r=  0  the  scheme  (35)  is  A-stable  in  the  sense  of  the  previous  deflnition, 
because  this  scheme  is  exact  for  any  exponential  function.  In  general,  the  characteristic  equation  (18) 


for  (34)  is 


P 


or 


1  + 


Xhr  P 


1+r 


1+r 


-  1 
1 


(P-1) 


Note  t  hat 


Xhr 

1+r* 


pr  * 1  iPi2(i+r)  -  PiPi5r  -  v\p? + y+t  -  pi+r + Pr  4-  p  - 1 

(/M)/+r- 1  *  (/+r  - 1  )(pi+r  - 1 ) 

Let  p  =»  Re,;?.  Then  the  real  part  of  the  numerator  is 
R2(l+r)  .  (Rl+Sr.R),.^  .  (R*+'-Rr)cosr0  -  1 


which,  for  r>0  and  R  >  1,  is  greater  than  or  equal  to 

R2(l+r)  .  (Rl+2r_R)  .  (RS+r.Rrj  .  j  _  (R.J)(R'.l)(R»+r.l)  >  0, 

whence  ReXh  >  0  in  (37).  Consequently,  if  r>0,  then  the  scheme  (34)  is  A-stable. 


(37) 


For  r<0,  ry^-1,  rewrite  (37)  as 

(1  -  £rXh )/>,+r  -  pr  -  p  +  (1  +  £rXh)  -  0.  '  (38) 

First,  assume  r  is  a  rational  number:  r  =*»  -  m/n,  where  m  and  n  are  two  positive  integers. 
Substituting  Z **pl/n  into  (38)  and  multiplying  by  Zm,  we  get 
Zn+m  -  (1  -  y^rXh)Zn  -  (1  +  y^rXh) Zm  +  1-0. 

Since  the  product  of  roots  of  the  above  equation  equals  (-])m+n,  if  any  root  is  greater  than  1  in 
magnitude,  then  at  least  one  root  is  less  than  1  in  magnitude.  By  continuity,  this  holds  for  all  real 
r<0,  including  r*-l.  Therefore,  if  r<0,  then  the  scheme  (34),  (35)  is  conditionally  A-stable. 


In  summary,  we  have  proven 


Theorem  21:  The  seheme  (84), (85)  it  A-stable  for  each  parameter  r  ^  0  and  it 

conditionally  A-ttable  for  r  <  0. 

Given  a  fixed  r,  scheme  (34)  is  exact  for  equations  y’  —  f(x,yfx))  «*=  Cj(x  -  C0)*/r,  where  CQ  and 
C,  are  constants,  just  as  the  trapezoidal  rule  is  exact  when  f  is  a  linear  function  of  x.  This  explains 
why  the  GMS  may  lead  to  better  results  near  a  singularity,  provided  we  can  find  a  good 
approximation  to  r. 

Note  that  we  can  interpret  the  GMS  as  an  Intergrand  Approximation  Method  (Jackson,  [4]). 
That  is,  the  discrete  numerical  solution  {YQ}  can  be  extended  to  a  continuous  approxmation  Y(x)  to 
the  solution  y(x)  of  (1)  satisfying  Y(xJ  y(xn)  by 

Y(x)  «  Y,  +  /  ‘  P(f(.,  Y(.));r](s)ds, 
where,  for  x  <  [x„.xn+1], 

Pbw  -  M\f 

h  a 

P[g:r]  is  a  nonlinear  interpolation  operator  computed  by  first  raising  g  to  the  power  r,  then 
performing  linear  interpolation,  and  finally  back  transforming  by  raising  the  interpolant  to  the  power 
1/r.  Of  course,  if  r~l,  P  reduces  to  a  linear  interpolation  operator. 

A  similar  technique  can  be  employed  to  enrich  a  piecewise  linear  space  to  solve  singular  two- 
point  boundary  value  problems,  using  either  the  finite  element  method  or  the  finite  difference 
method.  (See,  Jiachang  Sun  [10],  [11]). 


In  order  to  obtain  a  more  accurate  scheme,  we  set 


f  f 


h  . 


„2  'X*  X  +  ri 
s  ■  # 

From  (36),  the  scheme  (34)  is  fourth-order  accurate  for  this  value  of  r. 


(39) 


It  is  worth  mentioning  that  the  function  in  (39) 
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ff” 

F(f,rr )  -  pj-  (40) 

often  remains  bounded  even  when  f  and  its  derivatives  are  unbounded.  For  instance,  F  is  a  constant 
for  any  power  function  f»Cj(x  •  C0)\  And,  what  is  more  interesting,  F  is  identically  equal  to  1  for 
any  exponential  function,  f*Cj  Exp(C^x  -  CQ)). 

However,  the  F  is  not  easy  to  compute,  as  an  evaluation  of  f'  is  required.  Some  high-order 
schemes  based  on  non-polynomial  interpolation  developed  by  Lambert  [5]  and  Lambert  and  Shaw  [6] 
have  not  been  used  widely,  possibly  because  they  too  require  the  evaluation  of  higher  derivatives  of 
f.  Furthermore,  it  is  not  clear  that  these  methods  are  applicable  to  systems  of  equation. 


Fortunately,  we  can  avoid  computing  f 1  in  (39)  by  setting 


1  i  f"+> 

E<  — 

1  n+l 


L,. 

r 


(41) 


With  this  approximation,  the  scheme  (34)  retains  its  fourth-order  rate  of  convergence.  Also,  it 
retains  exact  for  f(x,y(x))  *  Cj(x  -  C0)*^r. 


Computing  experiments  show  that  there  is  only  a  slight  difference  in  accuracy  between  using 
(39)  and  (41)  in  the  scheme  (34);  sometimes  one  is  a  little  more  accurate,  and  sometimes  the  other. 
But  (41)  saves  computing  time,  and  f  is  not  required. 


An  alternative  derivation  of  the  GMS  is  obtained  by  taking 


Ff  -  f>A  (  r  /  0  ),  F  Jf  -  fr,  G(u)  -  /  *  F£  d£. 

®0 


With  this  notation,  (34)  may  be  written  as 


SMn+.J-r'fn+rF-'y  G(y) 


-If 


G^fn+ti-G^g 

F'*f  -  F'*f 
r  ‘n+l  r  'n 


(42) 


(43) 


This  formulation  can  be  generalized  by  considering  other  functions  F.  Using  Theorem  7  and  the 


local  truncatioD  error  formula  (11),  we  have 

Theorem  22:  For  any  one-to-one  map  F,  (48)  defines  a  second-order  one-step  implicit 
scheme  with  local  truncation  error 

L,f;  F)  _ .  _  { r  +  r*  +  0(1.*).  (44) 

This  formulation  unifies  most  of  schemes  described  in  this  and  the  previous  section.  Some 
examples  follow. 


1.  Ff 

2.  Ff 

3.  Ff 

4.  Ff 


f  leads  to  the  trapezoidal  rule. 


erf  leads  to  the  scheme  (35).  (It  is  independent  of  the  parameter  r). 


(l-f2)i/2 

af 

<T7fy/J 


leads  to  the  Elliptic  Scheme  (28). 


leads  to  the  Hyperbolic  Scheme  (30). 


Many  other  schemes  can  be  derived  using  this  formulation. 
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5.  Some  computational  considerations  for  the  GMS 

First,  we  consider  how  to  solve  the  implicit  equation  (2).  For  an  initial  value  system  (1)  with 
steep  gradients,  functional  iteration  may  be  employed: 


+  k  S(  f„,  ;  r,k  ,i ), 


where 


,»•!]  _  i  ,  f"+>  ^  ”  .  !l.  , 

<■  r. '' 

For  stiff  problems,  a  Newton  iteration  should  be  used  instead. 


As  a  simple  stopping  criterion  for  the  iteration,  we  use 

II  Y«t.  -  Y|k''U  II  < 

where  €y  is  a  parameter  to  be  specified. 

The  starting  value,  Y^n+j,  is  computed  by  a  conventional  explicit  method.  For  simplicity,  we 
use  the  Euler  Method  for  the  numerical  tests  in  the  last  section.  Of  course,  a  more  accurate 
predictor  may  be  used. 

The  rate  of  convergence  of  the  iteration  (45)  depends  upon  the  value  of  the  "contraction  factor” 

,ds  ,  asar  ^.asar 

h  ar*  h  Wd y  +  h  3F3y’  (47> 

where,  for  the  GMS, 

as  s  f'+w-f.'+W.  fW-rw, 


If  f  is  continuous,  then, 


S~>—  f„.  ash->0,f~>fn, 


§->  ^  +  +  K  ■  -  »• 

Also,  from  Theorem  5,  dS/df  -->  1/2  as  h  — >  0.  Therefore,  for  sufficiently  small  h,  iteration  (45)  is 


convergent  and,  moreover,  its  rate  of  convergence  is  close  to  that  of  the  trapezoidal  rule  using  the 
same  functional  iteration  procedure. 


It  is  worth  mentioning  that  the  GMS  is  particularly  well-suited  to  solving  problems  having  steep 
gradients,  especially  those  problems  for  which  f(x,y(x))  behaves  like  a  piecewise  power  function  of 
x.  In  practice,  it  may  be  more  efficient  to  use  the  GMS  only  on  the  sections  of  the  problem  having 
steep  gradients  and  a  conventional  scheme  on  the  sections  of  the  problem  where  f(x,y(x))  is  well- 
behaved.  We  consider  two  fourth-order  hybrid  schemes  of  this  type.  If  |fn|<£f*,  then  Mixl  uses  the 
cubic  Hermite  scheme  (modified  trapezoidal  rule) 


n+1 


Y  +  h  S±L 

n 


(fWU,)- 


(48) 


+  IT*1  n  *  '  n+l> 

* 

And  MixII  uses  the  classical  fourth-order  Runge-Kntta  scheme  in  place  of  (48).  If  |fj>f*,  then  both 
Mixl  and  MixII  use 


S<fu’fn+l) 


fn+.  -  fn 

Lo*  <f  A+.) 

f  -  r1’ 

*n+l  'n 


if  Irl^Ler  , 


if  |r+l|^.er*, 


and 


Slf  f  )  mm  r  fn+ A+l/f  Jr  '  fn 

”  “+l  m  <'„♦  A> '  •  >  ' 


(49) 


where  f*  and  er*  are  constants.  In  our  numerical  tests,  we  take  f*  — «  2  and  er*  »  0.01. 


6.  Extension  to  Systems 

So  far,  our  discussion  has  been  restricted  to  scalar  equations.  However,  we  have  used  nonlinear 
implicit  one-step  schemes  successfully  to  solve  systems  of  equations;  some  numerical  results  are 
presented  in  the  last  section. 

For  systems,  we  apply  the  scheme  (2)  to  the  individual  components  of  the  system.  The 
parameter  r  is  a  vector  whose  components  elements  are  determined  componentwise  by  formula  (46), 
where 

-  d  ,  u  dt  ,  dt 

r  -  -  37  +  E  ffir, 

is  a  vector.  Hence,  the  Jacobi  matrix  of  f(x)  is  needed  to  be  computed  for  finding  the  index  vector  r 
to  get  a  fourth  order  scheme.  The  advantage  of  (41)  over  (39)  is  more  significant  for  systems  than 
for  scalar  equations. 

The  analysis  of  nonlinear  schemes  for  systems  of  equations  is  an  open  problem  that  we  will 


consider  in  the  future. 
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7.  Numerical  Tests 

In  this  section,  we  present  some  numerical  results.  Throughout,  y(x)  denotes  the  solution  of  the 
problem  (1)  and  C(x-£K  is  a  piecewise  approximation  to  function  of  y’(x),  where  p  —  1/r  is  the  index 
of  the  approximation.  If  p  <  0,  the  approximate  position  of  the  singularity  is 

{-x.  +  MIp-l'-i)1. 

fn+l 

Throughout  this  section,  we  use  the  following  abbreviations: 

GMS  ~  the  scheme  (34); 

R-K  ~  the  classical  fourth-order  Runge-Kutta  scheme; 

C-H  —  cubic  Hermite  scheme  (48); 

Mixl  —  the  hybrid  scheme  composed  of  GMS  and  C-H  (49); 

MixII  —  the  hybrid  scheme  composed  of  GMS  and  R-K; 

L-S  --  the  scheme  proposed  by  Lambert  and  Shaw  (6); 

Error  —  Yn  -  y(x)  for  the  GMS; 

Er(f)  --  the  error  in  the  first  derivative  Y'  -  y’(x)  for  second-order  equations. 


A  uniform  mesh,  h  =  xn+J  •  xn,  is  used  throughout  the  section.  The  Fortran  program  was  run 
in  double- precision,  on  a  DEC-System  2060  computer  at  Yale.  The  iteration  error  cy  is  taken  to  be 
10'10. 


Test  1.  [6]  y’  =  1  +  y2,  y(0)  =  1 

with  solution  y(x)  =*  tan(x  +  ^)  which  has  a  strong  singularity  at  x 
Table  1.  x  =  0  (0.05)  0.75 


X  1 

y 

1  GMS 

1  P 

1 

L-S 

R-K 

0.701 

11.6814 

111.6808 

1  -1.975 

0.7828 

1 

11.6813 

11.6680 

0.751 

28.2383 

128.2305 

1  -1.992 

0.7851 

28.2378 

27.6947 

Remark:  For  the  exact  solution,  p  =  -2,  (  =  j  ~ 


0.7835. 


24 


Test  2.  (6]  xy’  —  y  +  5xV/Sx,  y(l)  —  0 

with  solution  y(x)  =  -5xLog(2-x)  which  has  a  weak  singularity  at  x  2. 
Table  2.  x  =  0  (0.05)  1.95 


X  1  y 

1 

GMS 

1  P 

Ban 

mm 

R-K 

1.901 

1 

21.8745 

21.8753 

1  -1.055 

2.001  1 
I. 

21.8748 

21.8746 

1.951 
.  1. 

29.2084 

29.2098 

1  -0.997 

2.001  1 
_ L 

29.2099 

29.2077 

Remark 

For 

the  exact 

solution. 

P  =  -1. 

e  =  2. 

Test  3.  [6]  (l-x)y’  ■■  y  Logy,  y(0)  —  e02 

with  solution  y(x)  —  e°-2/(M  which  has  an  essential  singularity  at  x  *»  1. 


Table  3.  x  =  0  (0.05)  0.95 


x  1 

y 

GMS  | 

P 

c 

1  L-S 

R-K 

6766 1  ~ 

7.3891 

7.3902  1 

-2.521 

0.963 

1  7.3954 

7.3646 

6755  T 

54.5982 

54.8956  1 

-3.126 

0.976 

1  57.1189 

47.1138 

Remark 

For 

the  exact  solution. 

ii 

i 

8 

*H 

ii 

*»/ 

These  results  show  that  the  GMS  is  more  accurate  than  the  classical  fourth-order  Runge-Kutta 


scheme  in  Tests  1  and  3;  the  accuracy  is  about  the  same  in  Test  2.  Also,  the  GMS  is  more  accurate 
than  the  scheme  of  Lambert  and  Show  in  Test  3;  the  accuracy  is  about  the  same  in  Tests  1  and  2. 
However,  the  GMS  doesn’t  need  f”’  which  the  L-S  scheme  requires. 
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Test  4.  (Artificial)  A  system  of  equations  consisting  of  four  components: 
y,’  —  -0.1  e'y2/y3,  y,(0)  —  .35 

>V  -  2{(y j-0.1 )  (lOy/}1/2,  y^O)  -  0 

y3’  *  e2y2  Log  y4,  y3(0)  —  0.1 

y4’  -  2y3y4,  y4(0)  -  e02. 

The  solution  is 

y,(x)  —  0.25(l-x)*  +  0.1,  y2(x)  —  -Log(l-x), 
y3(x)  «=  O.l(l-x)"2,  y4(x)  —  e®-2/(l-x)# 

Table  4  h  =  0.05 


X  1 

1 

yk(x) 

1  P 

1 

CHS 

R-K 

C-H 

0.9001 

(1)1 

0.1000 

1  3.01 

1 

0.77 

1 

1 

0.1000 

0.1000 

0.1000 

(2)1 

1 

2.3026 

1  -1.00 

1 

0.91 

1 

1 

2.3025 

2.1122 

2.2580 

(3)1 

1 

10.0000 

1  -3.00 

1 

0.97 

1 

1 

9.9992 

8.8429 

9.6617 

(4)  i 
_ 1 

7.3890 

1  -0.25 

1 _ 

0.96 

1 

1 

7.3867 

7.0125 

7.1183 

0.9501 

(1)1 

0.1000 

1-11.54 

1 _ 

0.27 

1 

1 

0.1000 

0.0999 

* 

(2)1 

1 

2.9957 

1  -0.50 

1 

0.95 

1 

1 

2.4771 

3.5606 

* 

(3)1 

1 

40.0000 

1  5.53 

1 . 

0.84 

1 

1 

36.4946 

40.0606 

* 

(4)1 

54.5982 
(  h/2 

1  -3.27 
=  0.025  ) 

0.98 

1 

1 

50.0038 

54.4607 

32.9154 

46.0766 

* 

44.9161 

*:  The  C-H  scheme  overflows  on  the  last  point  using  (45). 

Remark:  On  the  last  point  x=0.95  the  vector  of  first  derivatives  is  equal  to 

(  -0. 16x10"®,  2.1.  1 .4x10®,  3.8xl03  ). 


Test  5.  A  problem  with  an  integrable  singularity: 
v  1  ♦ 

y’  —  f+x  +  (i+*)|x  •  ,  y(o)  —  -  » 

The  solution  is 

y(x)  —  Sign(x  -  ^)  |x  -  j|I+t. 

Table  5.  Errors  t  -  0.1  (0.1)  1.0 


I  t=-0 . 1 


- 1  <  t  < 


Remark:  the  right  column  is  the  error  for  another  problem 

y’-d+tMx-O.S)*1,  y(0)-f^i,  -KK0, 
the  solution  of  which  is: 

y(x)  -  (0.5-x)1+t/(l+t),  for  x  <  1/2 

y(x)  *  (x-0.5)1+t/(l+t),  for  x  >  1/2 
which  has  a  turning  point  at  x  =  y. 

The  GMS  scheme  may  be  used  with  invariant  imbedding  to  solve  linear  two-point  boundary 
value  problems  with  various  singularity  properties.  For  these  problems,  y”  »  f  is  available  directly 
and  may  be  used  in  the  computation  of  the  fourth-order  GMS  scheme  using  (41)  to  compute  r. 

Consider  the  two-point  boundary  value  problem 

y”(x)  +  p(x)y’(x)  +  q(x)y(x)  —  r(x)  (50) 

aoy(°)  +  b^’fO)  —  c0,  a.y(l)  +  b,y’(l)  —  c.. 
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To  solve  this  problem,  we  use  the  sweep  method  of  Gelf&ad  and  Formin  [l]  (p.133).  (See,  also 
Miller  (7]  or  Scott  [8].) 


For  0. 

(i)  the  Initial  Value  Problem  for  the  Forward  Sweep  is 

u’  =  -  q  -  pu  -  u2,  u(0)  ■*  -  — , 

c0° 

v’  =  r  -  v(u  +  p),  v(0)  =  — . 

b„ 


(ii)  and  the  Initial  Value  Problem  for  the  Backward  Sweep  is 

vMi) 


u  y  +  v, 


y(i) 


aj+bj^l) 


For  b0  —  0, 

(i)  the  Initial  Value  Problem  for  the  Forward  Sweep  is 


u’«l  +  u(p  +  q«),  u(0)  =■  -  — , 

e 

v’-u(r  +  qv),  v(0)  —  — . 

*0 

(ii)  and  the  Initial  Value  Problem  for  the  Backward  Sweep  is 

CjU(l)-  b,v(l) 


uy 


y  +  v, 


y(i) 


a,u(l)  +  b, 


Test  6.  An  unstable  two  point  boundary  value  problem 
y”  -  165y’  -  2700y  +  4.95e15*  —  0, 
y’(0)  —  0.015,  y(l)  —  O.OOle15 
with  solution  y(x)  —  0.001  e15*. 
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Table  7.1 

Max irnum 

error  in 

Y  at  the  nodes 

e  1/h  1 

. . I__ 

R-K 

GMS 

Mi  xl 

C-H 

MixII 

0.1  1  10  1 

1 . L_ 

0.31+0 

-0.18-1 

-0.24-3 

0.79-3 

0.37-1 

1  20  1 

0.11-1 

0.75-4 

-0.12-4 

0.47-4 

0.99-3  f 

1  40  1 

0.45-3 

0.10-4 

-0.74-6 

0.29-5 

0.16-3 

1  80  1 

0.91-5 

0.23-5 

-0.46-7 

0.18-6 

0.33-4 

0.01  1  100  1 

0.21+0 

-0.24-2 

-0.19-4 

0.56-3 

-0.24-2  f 

1  200  1 

0.75-2 

0.65-5 

-0.12-5 

0.34-4 

0.10-4 

1  400  1 

0.30-3 

0.81-6 

-0.76-7 

0.21-5 

0.17-5 

1  800  1 

0.15-4 

0.19-6 

-0.47-8 

0.13-6 

0.34-6  ► 

0.00111000  1 

0.20+0 

-0.12-3 

-0.19-5 

0.54-3 

-0.12-3 

Table  7.2  Maximum  error  in  Y*  at  the  nodes 


0.38-3 


0.1  I  10  I 


0.34+1  -0.19+0 


0.88-2 


0.46+0 


I  20  I 
I  I 


0.12+0  0.83-3  -0.12-3 


0.52-3 


0.25-1 


I  400  I 
I  I 


I  800  I 
I  I 


0.00111000  I 


0.31-1  0.82-4  -0.76-5 


0.16-2  0.19-4  -0.48-6 


0.19+3  -0.1C+0  -0.19-2 


0.21-3 


0.13-4 


0.76-3 


0.29-3 


0.54+0  -0.12+0 
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tl 

»  ' 

Tab's  7.3 

CPU  TIME  in 

seconds  on 

a  DEC  20. 

f  1/h  1 

R-K 

GMS 

Mi  xl 

C-H  MixII 

V 

0.1  1  10  1 

0.08 

0.57 

0.48 

0.38  0.23 

1  20  1 

0.15 

0.52 

0.41 

0.38  0.21 

1  40  1 

0.28 

0.64 

0.54 

0.62  0.37 

1  80  1 

0.55 

1.15 

0.93 

1.12  0.66 

0.01  1  100  1 

1  83 

6.66 

5.63 

6.20  2.62 

L  • 

1  20?  1 

3  68 

7.49 

6.78 

7.83  4.69 

■  - ' 

1  400  1 

7.30 

12.98 

11.48 

13.96  7.92 

1  800  1 

14.63 

23.36 

21.76 

24.77  15.62 

21.82 

58.30 

47.74 

50 . 17__  23.40 

m 

Tabla  7 

ma 

• 

1  e  1 

mm 

GMS 

Mill 

C-H  MixII 

:  . 

10.1  1 

0.31*0 

-0.18-1 

-0.24-3 

0.79-3  0.37-1 

a 

Error  I  0.01  1 

0.21*0 

-0.24-2 

-0.19-4 

0.56-3  -0.24-2 

1  0.0011 

0.20*0 

-0.12-3 

-0.19-5 

0.54-3  -0.12-3 

1  0.1  1 

0.34*1 

-0.19*0 

0.38-3 

0.88-2  0.46*0 

• 

Er  (f)  1  0.01  1 

0.21*2 

-0.24*0 

-0.19-2 

0.57-1  -0.24*0 

1  0.0011 

0.19*3 

-0.12*0 

-0.19-2 

0.54*0  -0.12*0 

:■  • 
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